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Abstract 

In  many  areas  of  acoustical  imaging,  such  as  ultrasonic  non-destructive  evaluation  (NDE),  a 
realistic  calculation  of  ultrasonic  field  parameters  and  associated  elastic  wave  scattering 
requires  the  treatment  of  discontinuous,  layered  solids  in  complex  geometries.  These  facts 
suggest  the  need  for  an  accurate  and  geometrically  flexible  numerical  approach  for  the 
simulation  of  the  ultrasonic  field,  rather  than  reliance  on  semi-analytic  solutions. 

In  this  paper  we  present  an  approach  for  solving  the  elastic  wave  equation  in  discontinuous 
layered  materials  in  general  complex  geometries.  The  approach,  based  on  a  direct 
pseudospectral  solution  of  the  time-domain  elastodynamic  equations  consists  of  five  steps. 
The  first  step  decomposes  the  global  computational  domain  into  a  number  of  subdomains 
adding  the  required  geometrical  flexibility  to  the  method.  Moreover,  this  decomposition 
allows  for  efficient  parallel  computations,  hence  decreasing  the  computational  time.  The 
second  step  in  the  method  maps  every  subdomain  onto  the  unit  square  using  transfinite 
blending  functions.  With  this  curvilinear  mapping  the  elastodynamic  equations  can  be  solved 
to  spectral  accuracy,  and  furthermore,  complex  interfaces  can  be  approximated  smoothly, 
thus  avoiding  staircasing.  The  third  step  of  the  method  deals  with  the  evaluation  of  spatial 
derivatives  on  Chebyshev-Gauss-Lobatto  nodal  points  whithin  each  subdomain,  by  means  of 
a  pseudospectral  approach,  while  the  fourth  step  reconstruct  a  global  solution  from  the  local 
solutions  using  properties  of  the  equations  of  elastodynamics.  Each  subdomain  can  be 
prescribed  with  either  open,  physical  or  stress  free  boundary  conditions.  Boundary  conditions 
are  applied  by  means  of  characteristic  variables.  In  a  final  step,  the  global  solution  is 
advanced  in  time  using  a  fourth  order  Runge-Kutta  scheme.  Several  examples  of  elastic  wave 
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scattering  related  to  ultrasonic  NDE  are  presented  as  evidence  of  the  accuracy  and  flexibility 
of  the  proposed  computational  method. 
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I.  INTRODUCTION 


Numerical  solutions  of  the  elastic  wave  equation  are  needed  to  study  wave  propagation  in 
complex  distributions  of  material  for  which  analytical  solutions  do  not  exist.  Such  complex 
distributions  are,  in  particular,  found  in  industrial  materials  and  have  shown  to  be  a  challenge 
for  ultrasonic  non-destructive  evaluation  (NDE)  systems  in  which  the  main  interest  is  the 
recorded  transducer  signals  (A-  and  B-scans)  arising  due  to  elastic  scattering  from  material 
cracks®  Much  has  been  written  on  scattering  of  ultrasonic  waves  by  inhomogeneities®knd 
its  effect  on  ultrasonic  images®and  quantitative  NDE® 

Modeling  of  the  elastic  wave  equation  began  around  1970  with  the  finite  difference  (FD) 
method.  Alford  et  o/.^model  acoustic  scattering  with  a  higher  order  scheme  and  Kelly  et  a/.^ 
show  how  complex  interfaces  can  be  incorporated  into  the  simulations.  These  techniques 
handle  most  complex  material  geometries  but  are  limited  by  numerical  dispersion,  preventing 
the  modeling  of  waves  propagating  over  large  distances,  as  well  as  the  inability  to  impose  the 
correct  conditions  on  the  statevariables  across  material  interfaces.  In  the  early  eighties 


pseudospectral  (PS)  methods  were  introduced  to  enable  more  accurate  long  time  simulations 
of  acoustic  and  elastic  scattering.  Kosloff  et  o/.^and  Fornbcrg^solvc  the  acoustic  wave 


equation  using  a  two-dimensional  Fourier  PS  method  and  conclude  that  it  is  more  accurate 
than  a  FD  scheme.  Another  early  paper  of  interest  is  Ref.  [13  ,  which  shows  how  Fourier  PS 
methods  can  be  applied  to  problems  with  complex  interfaces.  The  Fourier  PS  method  is 
reviewed  in  Refs.  [  14  1&  .  Fourier  series  are  convenient  for  problems  with  periodic  boundary 
conditions,  but  when  the  solutions  are  non-periodic,  polynomial  approximations  are  a  more 
natural  choice.  Raggio1 — 'solves  the  acoustic  wave  equation  with  a  Chebyshev  PS  scheme. 
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Tessmer  et  o/.^and  Kosloff  et  a/.^apply  a  combined  Fourier  and  Chebyshev  method  to 
compute  wave  propagation  in  a  seismic  environment  with  surface  topography  and  propose  a 
three-dimensional  implementation^  Moreover,  essentially  one-dimensional  multidomain 
formulations  have  been  proposed  for  irregular  domain s^-^P  In  order  to  allow  elastic  waves  to 
pass  out  from  the  domain  without  reflections,  absorbing  boundary  conditions  for  elastic 
waves  have  been  used  in  many  variations.  See  for  example  Refs.  [23  27  .  Moreover, 

Rol 

mappings  may  be  used  to  enhance  the  accuracy  of  the  Chebyshev  PS  method  .  More 
recently,  Chebyshev  spectral  multidomain  techniques  have  become  a  standard  tool  in  fluid 
dynamics^0 J  and  is  emerging  as  such  in  computational  electromagnetic s A  more 
comprehensive  review  of  spectral  methods  for  hyperbolic  problems  can  be  found  in  Refs. 


In  this  paper,  we  present  an  approach  for  solving  the  elastic  wave  problem  in  general 
complex  distributions  of  materials  using  a  pseudospectral  multidomain  formulation.  The 
pseudospectral  elastodynamic  (PSE)  method  computes  a  direct  solution  to  the  elastodynamic 
equations  in  the  time  domain.  In  this  approach,  the  general  computational  domain  is  split  into 
a  number  of  smaller  subdomains,  each  chosen  such  that  they  can  be  smoothly  mapped  onto  a 
unit  square.  This  decomposition  is  performed  in  a  fully  bodyconforming  way  to  avoid 
problems  with  staircase  approximations  and  the  associated  errors.  This  enables  the 
representation  of  general  material  distributions  and  allows  for  the  construction  of  a 
pseudospectral  approximation  of  derivatives  within  each  of  the  smaller  domain  and  hence  an 
accurate  updating  of  the  local  fields.  To  recover  the  global  solution  from  the  many  local 
solutions  a  characteristic  decomposition  in  homogeneous  regions  of  the  computational 
domain  is  used  while  physical  boundary  and  interface  conditions  are  imposed  where  required. 


Filename:  EPS2D  rev3 


Page  5 


Nielsen,  JASA 


The  elastic  equations  are  advanced  in  time  using  a  4th  order  accurate  explicit  Runge-Kutta 
method  and  an  absorbing  matched  layer  is  used  to  truncate  the  computational  domain  in  an 
approximately  reflectionless  manner. 

In  the  following  section  the  governing  elastodynamic  equations  are  introduced  on  vector  form 
and  scattering  by  elastic  waves  are  discussed.  In  section  III,  the  features  of  the  multidomain 
Chebyshev  scheme  are  described.  Section  IV  gives  examples  of  wave  propagating  in  elastic 
half-spaces  illustrated  by  snapshots  of  the  velocity  field  at  particular  times  and  serves  as  an 
evaluation  of  the  scheme  while  Section  V  contains  a  few  concluding  remarks. 

II.  FORMULATION 


This  section  contains  two  parts.  In  part  A,  the  governing  elastodynamic  equations  are 
introduced  on  vector  form  and,  in  part  B,  elements  of  elastic  wave  scattering  are  reviewed. 

A.  Elastodynamic  Equations 

The  governing  elastodynamic  wave  equation  for  a  two-dimensional  isotropic  solid  are  based 
on  a  solution  of  the  equations  of  conservation  of  momentum  combined  with  the  stress-strain 

Eil 

relations  for  a  linear  elastic  solid  undergoing  infinitesimal  deformations  .The  elastodynamic 

CT71 

equations  are  given  by  a  system  of  two  coupled  wave  equations,  as 
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In  this  system  of  equations  u=(ux,  uz)T  is  the  displacement  vector,  o=(oxx,  azz,  axz )  represents 
the  symmetric  stress  tensor  and  (x,z)  are  the  Cartesian  coordinates.  A(x,z)  and  ju(x,z)  are  the 
Lame  constants  (i.e.  rigidity  and  shear  modulus,  respectively,  for  fluids:  ju=0),  t  is  the  time 
and  p(x,z)  is  the  mass  density.  The  body  force  is  given  as  f=(fx,fz)  and  represents  the  applied 
source. 

The  elastodynamic  equations  in  (1)  can  be  recast  into  a  hyperbolic  velocity- stress  system  of 
first  order  in  time.  This  formulation  consists  of  five  coupled  first-order  partial  differential 
equations 
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Here,  v  is  the  particle  velocity  vector  related  to  u  as  v  =  — . 

dt 

Introducing  q  =  [vv,vz,ara,(T__,a'tz]r  as  the  state  vector,  describing  the  state  of  the  system,  Eq. 
(2)  takes  the  simple  vector  form 


?1  =  A?1+A  ^L+s 

dt  ldx  2  dz 


(3) 


where  A  i  and  Ai  are  matrices  containing  the  isotropic  material  parameters 
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0  , 

and  the  body  forces  are  given  as  S  =  [~pfx, ~pfz, 0,0, 0]r  . 


B.  Scattering  by  Elastic  Waves 


In  general,  scattering  originates  from  waves  of  different  type  impinging  on  complex 
interfaces.  The  most  significant  waves  are  compressional  waves  and  vertically  polarized 
shear  waves  (denoted  P-waves  and  S-waves,  respectively)  with  wave  speeds 
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c„ 


a +2  ju) 


The  ratio  between  the  velocities  of  the  two  wave  types  is 


_Q  =  I  l-2v 
Cp  \  2(1- v)  ’ 


(5) 


(6) 


where  the  Poisson  ratio,  v,  being  given  as 


1 

v  =  — 
2 


1-2  C)IC] 

i-c;/c2p 


(7) 


For  industrial  solids  related  to  ultrasonic  NDE,  it  is  often  assumed  that  the  Poisson  ratio  is 
between  0-0.5.  In  ferritic  steel,  for  example,  the  Poisson  ratio  is  around  0.29  and  thus 
Cs  /  Cp~  0.5  .  This  means,  that  the  compressional  velocity  is  about  50  percent  larger  than  the 

shear  wave  velocity.  Since  the  PSE  approach  is  based  on  the  elastodynamic  wave  equations 
without  physical  approximations,  the  approach  accounts  not  only  for  direct  P-  and  S-waves, 
reflected  waves  (PP,  SS),  and  multiply  reflected  waves  (PPP,  SSS),  but  also  for  converted 
reflected  waves  (PS,  SP),  refracted  waves  (PR,  SR),  diffracted  waves  (PD,  SD),  head  waves 
(H),  and  interfacial  waves.  Among  the  latter  are  for  example  Rayleigh  (R)  and  Stoneley  (St) 
waves.  Rayleigh  waves  travel  along  the  surface  of  the  solid  while  Stoneley  waves  travel 
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along  internal  interfaces  between  two  elastic  domains.  Generally,  these  waves  travel  slower 
than  the  shear  wave  and  decays  exponentially  away  from  the  interface. 

III.  NUMERICAL  APPROACH 

The  PSE-approach  for  solving  the  elastic  wave  equation  in  complex  geometries  containing 
general  material  properties  consists  of  a  few  steps.  The  first  step  decomposes  the  global 
computational  domain  into  a  number  of  bodyconforming  subdomains.  The  second  step  in  the 
method  maps  every  subdomain  onto  the  unit  square  using  transfinite  blending  functions. 
These  two  steps  enable  the  elastodynamic  equations  to  be  evaluated  using  pseudospectral 
methods,  and  furthermore,  complex  interfaces  can  be  approximated  in  a  smooth 
bodyconforming  way.  The  third  step  deals  with  the  construction  of  a  global  solution  from  the 
local  solutions  as  well  as  enforcing  the  correct  boundary  condition  as  being  either  open, 
physical  or  stress  free.  In  the  final  step,  the  global  solution  is  advanced  in  time  using  a  fourth 
order  Runge-Kutta  scheme. 

This  section  contains  a  discussion  of  the  PSE-approach.  In  part  A,  complex  geometries  are 
incorporated  using  a  curvilinear  representation.  In  part  B,  the  approximation  of  spatial 
derivatives  by  a  Chebyshev  collocation  scheme  is  discussed,  while  part  C  addresses  how 
local  solutions  are  patched  using  characteristic  variables  to  recover  the  global  solution. 
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A.  Curvilinear  Formulation 

To  enable  the  representation  of  general  distributions  of  materials  and  interfaces  while 
maintaining  a  high  accuracy,  the  global  computational  domain  is  split  into  a  number  of 
smaller  subdomains,  each  chosen  such  that  they  can  be  smoothly  mapped  to  a  unit  square.  As 
the  general  blocks  are  curvilinear  a  mapping  is  chosen  to  connect  the  physical  grid  with  the 
local  computational  grid  as  illustrated  in  Fig.  1.  The  curviliniar  physical  grid  has  the 
coordinates  (x,z)  whereas  the  rectangular  grid  has  the  coordinates  (£  ?])  connected  as 

g  =  g(x,z),  J)  =  tj(x,z).  (8) 


By  applying  the  chain  rule,  Eq.  (3)  can  be  written  in  the  curvilinear  representation  as 


^-  =  A(V£)||  +  A(V/7)|*-  +  S.  (9) 

dt  dg  dr/ 


Setting  n=(nx,nz),  the  matrix  A(n)  is  given  as 


0  0  nj  p  0  njp " 

0  0  0  njp  njp 

A(n)  =  (A +  2  p)nx  An.  0  0  0.  (10) 

Anx  (A  +  2p)n.  0  0  0 

v  0  0  0  , 
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Transfinite  blending  functions  are  used  to  establish  a  connection  between  the  physical 
curvilinear  grid  and  the  auxiliary  gird.  A  complete  treatment  of  these  functions  may  be  found 
in  Gordon  et  a/.^8^9J 

B.  Pseudospectral  Differentiation 

In  the  Chebyshev  method  a  function  q(£)  is  approximated  by  the  polynomial 

7=1 


that  interpolates  q(c,)  at  N  distinct  collocation  points  £ . .  Due  to  the  non-periodic  nature  of  the 
problem  the  grid  points  are  chosen  as  the  Chebyshev-Gauss-Lobatto  points 


€j  =  cos 


l  N- 1  J 


7  =  1,.~,  N. 


(12) 


and  interpolation  between  the  collocation  points  is  based  on  the  Chebyshev-Lagrange 
polynomials  given  as 


CjiN-lfg-Zj) 


(13) 


where  c,  =  cN  =  2  and  c2  =  c3  = ...  =  cN_  (  =  1 .  Tn.j  is  the  Chebyshev  polynomial  of  order  N-l 
defined  as 
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TN-i (£)  =  cos((N - 1) COS-1  £,)  ,  (14) 

where  the  computational  domain  is  the  region  <f  e  [-1,1]. 

The  spatial  derivative  of  q(Q  is  calculated  by  analytically  differentiating  the  interpolation 
polynomial  and  the  differentiation  operator  may  then  be  represented  by  a  matrix 
Djk  =  </>k(J;.).  Hence,  the  approximative  derivative  of  q(Q  at  the  collocation  points  can  be 
found  by  a  matrix  D  multiplying  the  values  of  q(Q  at  the  collocation  points 

as) 

k= 1 


where  the  differentiation  matrix  has  the  entries29,34 


D 


jk 


ck  (-1  )1+k 

Cj 

1  & 

2  (1-4  )2 

2(V-lf  +1 
6 

2(iV-l)2+l 

6 


j*k 

j  =  k  ^l,N 

j  =  k  =  1 
j  =  k  =  N 


(16) 


2 

The  cost  of  computing  the  derivative  from  the  matrix  multiplication  is  0(/V~ )  operations.  In 
two  dimensions,  spatial  derivatives  are  calculated  using  Eq.  (18)  and  a  matrix-matrix  product. 
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C.  Multidomain  Formulation  and  boundary  conditions 

Once  the  computational  domain  is  divided  into  a  number  of  subdomains  consisting  of 
curvilinear  quadrilaterals,  the  spatial  derivatives  are  computed  locally  on  each  domain,  and 
the  solution  to  the  elastodynamic  equations  are  advanced  in  time  with  a  fourth-order  Runge- 
Kutta  scheme.  In  order  to  recover  a  global  solution  from  the  subdomains,  two  strategies  of 
interchanging  internal  information  are  employed.  On  boundaries  between  materials  with 
different  elastic  constants,  continuity  of  stress  and  displacement  vector  is  required.  On 
boundaries  between  materials  with  the  same  elastic  constants  a  reconstruction  based 
characteristic  variable  is  employed.  The  boundary  conditions  at  a  stress-free  surface  require 
the  stress  to  vanish. 

In  order  to  patch  the  local  solutions,  the  characteristic  variables  of  the  matrix  A  =  Ajnx  + 
A 2fiz  =  QAQ1  is  used.  Hence,  the  diagonal  eigenvalue  matrix  An=X,  has  the  entries 

1  /V  /V 

corresponding  to  the  characteristic  velocities  (0,Cp,-Cp,Cs,-Cs)  and  Q  with  n=\n\(nx,nz )  is 
given  as 


0 

0 

-2  jun^-A(n^-n^) 

A{A+/i)nzn} 

A+2  fi 

A+2/u 

A+2  fi 

nx\j  A+2  fi 

nz-\J  A+2fi 
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nz 
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2  nxnz 
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nx :\j  A+2p 
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A+2  n 
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A+2/u 

2  nxnz 
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—\nzy[p/4 

2Ux-JPM 

~\nxnz 

\nxn. 

1  (»;-»; 

~2nx\[PM 

~\nxnx 

\nxnz 

I  (»;-»; 

(17) 


The  characteristic  variables  take  the  form 
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R  =  Q-\n)q  =  [Ri  R2  R3  R4  Rsf .  (18) 

Once  the  characteristic  variables  are  recovered,  one  can  exploit  the  signs  of  the  corresponding 
eigenvalues  to  determine  how  to  pass  information  between  two  neighboring  elements. 
Indeed,  based  on  the  signs  of  the  eigenvalues,  one  may  observe  that  R2  and  R4  will  be 
propagating  along  the  direction  of  the  normal,  i.e.,  they  will  leave  the  computational  domain. 
Hence,  they  need  not  to  be  altered  but  will  enter  the  neighboring  element.  Conversely,  R3  and 
R5  will  enter  the  current  element  and  will  need  to  be  prescribed.  This  is  done  by  passing  the 
information  from  R2  and  R4  from  the  neighboring  element  into  the  present  one  as  what  leaves 
one  domain  must  enter  the  neighboring  domain.  Finally  it  may  be  observed  that  R\  is  not 
propagating  and  it  is  required  to  be  continuous.  Once  the  characteristic  variables  have  been 
updated  within  the  element  as  discussed  above,  the  corrected  state  vector  can  be  recovered 
and  the  patching  of  two  elements  is  completed.  This  procedure  is  then  repeated  along  every 
element  interface  at  every  time  step  to  continuously  recover  the  global  solution  from  the 
many  local  solutions  and  ensure  that  information  is  propagating  in  accordance  with  the  nature 
of  the  problem. 

Open  boundary  conditions  for  simulating  an  infinite  medium  and  avoiding  non-physical 
reflections  from  the  boundaries  of  the  numerical  grid  is  implemented  using  a  matched 
absorbing  layer.  This  layer  is  introduced  through  an  absorbing  term  added  to  Eq.  (9)  as 
described  in  Ref.  [  W  . 
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IV.  NUMERICAL  RESULTS 


In  this  section  numerical  examples  of  elastic  wave  scattering  related  to  ultrasonic  NDE  are 
presented.  In  part  A,  elastic  scattering  from  a  plane  slit  (crack)  is  calculated.  Part  B  discusses 
elastic  scattering  from  a  two-layer  elastic  solid,  while  part  C  presents  snapshots  of  scattering 
by  an  elastic  cylinder. 

A.  Elastic  Scattering  from  a  Slit 

The  first  example  involves  wave  propagation  in  a  two-dimensional  homogeneous  elastic  half¬ 
space  with  a  body  force  applied  at  the  stress-free  surface.  This  example  is  a  classical 
modeling  problem  in  NDE  called  Lamb’s  problem.  The  P-  and  S-wave  velocities  for  the  solid 

are  set  to  1  and  V3,  respectively,  corresponding  to  a  Poisson  ratio  of  0.25.  As  a  first  test 
case,  the  2D  pseudospectral  code  is  compared  with  a  3D  analytical  solution.  The  analytical 
solution  was  found  by  the  Cagniard-de  Hoop  formalism  giving  an  analytical  solution  in  time 
for  Green’s  function.  The  corresponding  velocity  component  was  then  derived  by 
convoluting  the  Green’s  function  with  a  source  signal,  which  here  was  set  to  a  pulsed  raised- 
cosine  signal.  The  numerical  and  analytical  time  solutions  (A-scans)  were  compared  for  vz. 
From  Table  I  it  can  be  seen  that  there  is  excellent  agreement  between  the  two  solutions  even 
for  as  few  as  14  grid  points  per  spatial  wavelength.  The  global  error  (or  infinity  norm)  was 
not  improved  for  decreasing  time  step,  due  to  the  fundamental  difference  between  the  2D  and 
3D  solution. 
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In  the  next  test  case,  the  body  force  is  a  directional  force  applied  at  the  stress  free  surface  and 
representing  a  source  transducer.  In  this  example  a  Gaussian  distribution  for  the  point  source 
was  used  and  a  time  history,  s(t)  set  to 

s(t)  =  -acos{lK{t-t0))e~^-to)2 ,  (19) 

where  a,  (5,  and  t0  are  constants.  The  current  examples  consider  a  rectangular  computational 
domain  consisting  of  40  subdomains,  each  with  16x16  grid  points,  as  shown  in  Fig.  2.  The 
source  was  set  on  a  stress  free  boundary,  while  the  other  sides  were  open  (absorbing) 
boundaries.  The  vertical  slit  illustrates  an  infinitesimal  thin  stress  free  notch  or  crack. 

Figure  3  displays  elastic  scattering  from  a  vertical  stress  free  slit  due  to  the  vertical  point 
force.  The  vertical  velocity  component  vz  is  seen  at  different  times  (a)  t  =  1.8  (is,  (b)  t  =  4.1 
(is,  and  (c)  t  =  5.5  (is.  Figure  3a,  shows  the  P-,  S-  and  R-wave  just  after  leaving  the  point 
source  at  the  free  surface  and  Figs.  3b  and  3c  illustrate  the  Head  (H)  wave  propagation.  The 
head  wave  originates  at  the  vertical  slit  together  with  the  P-wave  and  ends  at  the  S-wave  front 
where  it  is  perpendicular  with  the  radius  vector.  The  angle  of  the  head  wave  with  respect  to 
vertical  depends,  as  expected,  on  the  ratio  between  the  two  velocities,  i.e.  the  Poisson  ratio,  as 
seen  in  Eq.  (7).  Notice,  that  no  artificial  reflections  are  observed  at  the  open  boundaries  due 
to  the  matched  absorbing  layer. 

Figure  4  displays  elastic  scattering  from  a  horizontal  slit  due  to  the  vertical  point  force.  The 
vertical  velocity  component  v-  is  seen  at  different  times  (a)  t  =  4.1  ps,  (b)  t  =  5.5  (is,  and  (c)  t 
=  8.2  ps.  Figure  4a  shows  the  reflected  PP-wave  after  impinging  on  the  horizontal  stress-free 
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slit.  In  Fig.  4b,  the  wave  develops  and  several  wave  phases  are  present,  including  the 
diffracted  pressure  wave  (PD)  and  Rayleigh  wave  (R).  As  can  be  appreciated,  the  P-wave 
looses  amplitude  with  time  due  to  geometrical  spreading,  while  the  Rayleigh  wave  keeps  the 
same  amplitude  since  it  is  confined  to  the  slit.  Figure  4c  shows  the  multiply  reflected  wave 
(PPP)  due  to  the  stress-free  slit  and  the  free  surface. 

The  corresponding  time  history  (A-scan)  is  seen  in  Fig.  5  for  two  different  receiver  positions 
on  the  stress-free  surface.  The  first  receiver  position  is  identical  with  the  transmitter  (i.e.  a 
pulse-echo  transducer)  and  indicated  with  a  dotted  line.  In  this  position,  the  initial  wave  pulse 
is  followed  by  the  PP-wave  reflected  from  the  horizontal  slit.  The  second  receiver  position  is 
illustrated  in  Fig.  2  and  as  indicated  with  a  solid  line  in  Fig.  5,  the  S-wave  is  proceeded  by  the 
PP-wave. 

B.  Scattering  from  a  Two-Layer  Elastic  Solid 

The  second  example  illustrates  scattering  from  a  two-layered  elastic  solid.  The  model 
consists  of  two  different  elastic  regions,  e.g.  a  steel  rod  embedded  in  a  lead  coating  with 
material  parameters  given  in  Table  II.  The  calculations  were  performed  on  the  same  grid  as 
used  in  the  previous  example,  but  with  the  extension  of  an  elastic  steel  layer  in  5  subdomains. 
The  time  interval  was  19.5  ns  and  the  total  calculation  time  was  9.8  ps. 

Figure  6  illustrates  snapshots  at  times  (a)  t  =  7.4  ps  and  (b)  t  =  9.8  ps  of  the  particle  velocity. 
It  may  be  noticed  from  Fig.  6a,  that  the  P-wave  is  refracted  according  to  Snell’s  law  within 
the  steel  rod  and  that  the  polarity  of  the  reflected  P-wave  (PP)  is  changed.  In  Fig.  6b,  a 
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diffracted  P-wave  (PD)  is  seen  near  the  edge  of  the  steel  rod  and  a  Stoneley  (St)  wave  may 
also  be  identified  along  the  interface  of  the  two  elastic  solids. 

C.  Scattering  from  an  Elastic  Cylinder 

The  final  example  illustrates  scattering  by  an  elastic  cylinder.  This  example  is  chosen  because 
it  appears  frequently  in  the  NDE  community.  The  model  consists  of  the  same  two  different 
elastic  regions  as  used  in  the  previous  example.  The  transmitting  line  source  (i.e.  transducer) 
is  located  at  the  stress  free  surface  and  consists  of  1 1  point  sources  distributed  evenly  over 
three  subdomains.  The  computational  domain  consisted  of  29  subdomains,  each  with  20  x  20 
grid  points,  as  shown  in  Fig.  7.  Figure  8  displays  elastic  scattering  by  an  elastic  cylinder  due 
to  the  transducer.  The  velocity  component  is  seen  at  different  times  (a)  t  =  4.9  ps  and  (b)  t  = 
7.1  ps.  Figure  8a,  shows  the  P-wave  impacting  the  elastic  steel  cylinder.  The  steel  cylinder 
has  a  magnifying  effect  on  the  P-wave  due  to  the  faster  velocity  compared  with  the  lower 
velocity  in  lead.  In  Fig.  8b,  the  Stoneley  (St)  wave  is  seen  clearly  propagating  along  the 
cylinder.  For  this  calculation,  the  time  interval  was  2.37  ns  and  the  wave  propagation  time 
was  14.2  ps.  The  CPU-time  for  each  time  step  was  about  5.1  seconds,  corresponding  to  a 
total  computational  time  of  about  8  hours.  The  corresponding  backscattered  velocity  profile 
(B-scan)  was  calculated  from  11-point  receivers  distributed  over  the  aperture  of  the 
transducer.  The  B-scan  is  seen  in  Fig.  9  and  is  usually  the  main  component  in  solving  the 
inverse  problem.  Here,  it  corresponds  to  a  single  projection. 
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V.  CONCLUSIONS 

A  multidomain  pseudospectral  method  is  proposed  for  the  simulation  of  scattering  from 
complex  interfaces  in  elastic  media.  The  method  is  based  on  a  Chebyshev  pseudospectral 
formulation  and  computes  a  direct  time-domain  solution  of  the  elastodynamic  equations. 
Complex  interfaces  are  resolved  with  a  curvilinear  multidomain  representation  and  local 
solutions  are  patched  using  characteristic  variables. 

The  first  example  demonstrate  the  ability  of  the  PSE-approach  to  accurately  compute 
scattering  from  a  plane  stress  free  slit  and  illustrate  different  wave  types.  The  second  and 
third  example  demonstrate  elastic  scattering  from  two  elastic  regions  with  a  plane  and 
cylindrical  shape.  Clear  indications  of  interfacial  waves  caused  by  the  elastic  cylinder  were 
observed. 
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TABLE  I.  Global  error  of  the  velocity  component  y,  for 
Lamb’s  problem  as  a  function  the  grid  points  per  wavelength. 


Parameter 

Symbol 

Grid  Points  (N) 

12  x  12 

16  x  16 

16  x  16 

Grid  points 

Nppw 

10.6 

14.2 

14.2 

AV  avelength 

Step  size 

dt 

0.033 

0.018 

0.0078 

Global  error 

E 

0.13 

0.077 

0.077 
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TABLE  II.  Material  parameters  used  in  the  calculation  of  two-layer  simulation. 


Parameter 

Symbol 

Steel 

Value 

Lead 

P-wave  velocity 

Cp  (m/s) 

5900 

2050 

S-wave  velocity 

Cy  (m/s) 

3200 

670 

Poisson  ratio 

V 

0.29 

0.44 

Lame  constant 

X  (Gpa) 

113.2 

16.5 

Lame  Constant 

P  (Gpa) 

80.9 

5.5 

Mass  density 

p  (g/cm3) 

7.9 

11.3 
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Figure  Captions 


FIG.  1.  Mapping  between  physically  curved  grid  onto  auxiliary  rectangular  grid. 

FIG.  2.  Decomposition  of  Lamb’s  problem  into  40  subdomains. 

FIG.  3.  Elastic  scattering  from  a  vertical  slit  due  to  a  point  force  (a)  t  =  1.8  ps,  (b)  t  =  4.1  ps,  and 
(c)  t  =  5.5  ps. 

FIG.  4.  Elastic  scattering  from  a  horizontal  slit  due  to  a  point  force  (a)  t  =  4.1  ps,  (b)  t  =  5.5  ps,  and 
(c)  t  =  8.2  ps. 

FIG.  5.  Time  history  (A-scan)  comparison  between  elastic  v,  component  at  (.r,z)=(0,0)  and  (0,2). 
FIG.  6.  Snapshot  of  particle  velocity  in  two  different  elastic  solids  (a)  t  =  7.4  ps,  and  (b)  t  =  9.8  ps. 
FIG.  7.  Decomposition  of  cylinder  into  29  subdomains. 

FIG.  8.  Scattering  from  elastic  cylinder  due  to  a  line  source  (a)  t  =  4.9  ps,  and  (b)  t  =  7.1  ps. 

FIG.  9.  B-scan  of  backscattered  velocity  component. 
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FIG.  1.  Mapping  between  physically  curved  grid  onto  auxiliary  rectangular  grid. 
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FIG.  6.  Snapshot  of  particle  velocity  in  two  different  elastic  solids, 
(a)  t  =  1 A  (is,  and  (b)  t  =  9.8  (as. 
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Open  boundary 


FIG.  7.  Decomposition  of  cylinder  into  29  subdomains. 
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Transducer 


(a) 


FIG.  8.  Scattering  from  elastic  cylinder  due  to  a  line  source, 
(a)  t  =  4.9  (Us,  and  (b)  t  =  7.1  jus. 
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FIG.  9.  B-scan  of  backscattered  velocity  component. 
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